package net.wesc.ac.server;

public class LMPP implements LMfunc
{
  static final int a0 = 0; //beta0
  static final int a1 = 1; //beta1
  static final int a2 = 2; //C psi
  static final double b = 1.2;
  static final double alpha = 2.0;
  static final double A0 = 0.509;
  static final double B0 = 3.291;
  static final double temp0 = 298.15;
  static int Zc;
  static int Za;
  static int Va;
  static int Vc;
  static int V;
  static double temp = temp0;
  static double A = A0;
  static double B = B0;
  
  public LMPP(int tZc, int tZa, double ttemp){
	  Zc = tZc;
	  Za = tZa;
	  Vc = Za;
	  Va = Zc;
	  V = Va+Vc;
	  temp = ttemp;
	  A = A0*Math.pow(temp0/temp, 1.5);
	  B = B0*Math.pow(temp0/temp, 0.5);
  }
  
  public double[] initial()
  {
    double[] a = new double[3];
    a[a0] = 1;
    a[a1] = 1;
    a[a2] = 1;
    return a;
  } //initial

  public double val(double[] x, double[] a){
	  double I = 0.5*(Va*x[0]*Math.pow(Za, 2)+Vc*x[0]*Math.pow(Zc, 2)); //I ionic strength
	  double fgamma = -(A/3)*(Math.sqrt(I)/(1+b*Math.sqrt(I))+(2/b)*Math.log(1+b*Math.sqrt(I)));
      double Bgamma = 2*a[a0]+((2*a[a1])/(Math.pow(alpha, 2)*I))*(1-(1+alpha*Math.sqrt(I)-0.5*Math.pow(alpha, 2)*I)*Math.exp(-alpha*Math.sqrt(I)));
	  return Math.abs(Zc*Za)*fgamma + x[0]*(2*Va*Vc/V)*Bgamma + Math.pow(x[0], 2)*(2*Math.pow(Va*Vc, 1.5)/V)*a[a2];
  } //val
  
  public double grad(double[] x, double[] a, int a_k){
	  if (a_k == a0){
		  return 2*(2*Va*Vc/V)*x[0];
	  }
	  if (a_k == a1){
		  double I = 0.5*(Va*x[0]*Math.pow(Za, 2)+Vc*x[0]*Math.pow(Zc, 2)); //I ionic strength
		  return x[0]*(2*Va*Vc/V)*(2*Va*Vc/V)*(2/(Math.pow(alpha, 2)*I))*(1-(1+alpha*Math.sqrt(I)-0.5*Math.pow(alpha, 2)*I)*Math.exp(-alpha*Math.sqrt(I)));
	  }
	  if (a_k == a2){
		  return Math.pow(x[0], 2)*(2*Math.pow(Va*Vc, 1.5)/V);
	  }
	  else {
		  assert false;
		  return 0.;
	  }
  } //grad
} //DHTest
